home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / newmat.exe / NEWMAT2.CXX < prev    next >
C/C++ Source or Header  |  1991-07-30  |  11KB  |  343 lines
  1. //$$ newmat2.cxx      Matrix row and column operations
  2.  
  3. // Copyright (C) 1991: R B Davies and DSIR
  4.  
  5.  
  6. #include "include.hxx"
  7.  
  8. #include "boolean.hxx"
  9.  
  10. typedef double real;                 // all floating point double
  11.  
  12. #include "newmat.hxx"
  13.  
  14. #include "newmatrc.hxx"
  15.  
  16. //#define REPORT { static ExeCounter ExeCount(__LINE__,2); ExeCount++; }
  17.  
  18. #define REPORT {}
  19.  
  20. //#define MONITOR(what,storage,store) \
  21. //   { cout << what << " " << storage << " at " << (long)store << "\n"; }
  22.  
  23. #define MONITOR(what,store,storage) {}
  24.  
  25. /************************** Matrix Row/Col functions ************************/
  26.  
  27. void MatrixRowCol::Add(const MatrixRowCol& mrc)
  28. {
  29.    REPORT
  30.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  31.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  32.    if (l<=0) return;
  33.    real* elx=store+f; real* el=mrc.store+f;
  34.    while (l--) *elx++ += *el++;
  35. }
  36.  
  37. void MatrixRowCol::AddScaled(const MatrixRowCol& mrc, real x)
  38. {
  39.    REPORT
  40.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  41.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  42.    if (l<=0) return;
  43.    real* elx=store+f; real* el=mrc.store+f;
  44.    while (l--) *elx++ += *el++ * x;
  45. }
  46.  
  47. void MatrixRowCol::Sub(const MatrixRowCol& mrc)
  48. {
  49.    REPORT
  50.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  51.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  52.    if (l<=0) return;
  53.    real* elx=store+f; real* el=mrc.store+f;
  54.    while (l--) *elx++ -= *el++;
  55. }
  56.  
  57. void MatrixRowCol::Inject(const MatrixRowCol& mrc)
  58. // copy stored elements only
  59. {
  60.    REPORT
  61.    int f = mrc.skip; int l = f + mrc.storage; int lx = skip + storage;
  62.    if (f < skip) f = skip; if (l > lx) l = lx; l -= f;
  63.    if (l<=0) return;
  64.    real* elx = store+f; real* ely = mrc.store+f;
  65.    while (l--) *elx++ = *ely++;
  66. }
  67.  
  68. real DotProd(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  69. {
  70.    REPORT                                         // not accessed
  71.    int f = mrc1.skip; int f2 = mrc2.skip;
  72.    int l = f + mrc1.storage; int l2 = f2 + mrc2.storage;
  73.    if (f < f2) f = f2; if (l > l2) l = l2; l -= f;
  74.    if (l<=0) return 0.0;
  75.    
  76.    real* el1=mrc1.store+f; real* el2=mrc2.store+f;
  77.    real sum = 0.0;
  78.    while (l--) sum += *el1++ * *el2++;
  79.    return sum;
  80. }
  81.  
  82. void MatrixRowCol::Add(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  83. {
  84.    int f = skip; int l = skip + storage;
  85.    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
  86.    if (f1<f) f1=f; if (l1>l) l1=l;
  87.    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
  88.    if (f2<f) f2=f; if (l2>l) l2=l;
  89.    real* el = store + f;
  90.    real* el1 = mrc1.store+f1; real* el2 = mrc2.store+f2;
  91.    if (f1<f2)
  92.    {
  93.       register int i = f1-f; while (i--) *el++ = 0.0;
  94.       if (l1<=f2)                              // disjoint
  95.       {
  96.      REPORT                                // not accessed
  97.          i = l1-f1;     while (i--) *el++ = *el1++;
  98.          i = f2-l1;     while (i--) *el++ = 0.0;
  99.          i = l2-f2;     while (i--) *el++ = *el2++;
  100.          i = l-l2;      while (i--) *el++ = 0.0;
  101.       }
  102.       else
  103.       {
  104.          i = f2-f1;    while (i--) *el++ = *el1++;
  105.          if (l1<=l2)
  106.          {
  107.             REPORT
  108.             i = l1-f2; while (i--) *el++ = *el1++ + *el2++;
  109.             i = l2-l1; while (i--) *el++ = *el2++;
  110.             i = l-l2;  while (i--) *el++ = 0.0;
  111.          }
  112.          else
  113.          {
  114.             REPORT
  115.             i = l2-f2; while (i--) *el++ = *el1++ + *el2++;
  116.             i = l1-l2; while (i--) *el++ = *el1++;
  117.             i = l-l1;  while (i--) *el++ = 0.0;
  118.          }
  119.       }
  120.    }
  121.    else
  122.    {
  123.       register int i = f2-f; while (i--) *el++ = 0.0;
  124.       if (l2<=f1)                              // disjoint
  125.       {
  126.      REPORT                                // not accessed
  127.          i = l2-f2;     while (i--) *el++ = *el2++;
  128.          i = f1-l2;     while (i--) *el++ = 0.0;
  129.      i = l1-f1;     while (i--) *el++ = *el1++;
  130.      i = l-l1;      while (i--) *el++ = 0.0;
  131.       }
  132.       else
  133.       {
  134.          i = f1-f2;    while (i--) *el++ = *el2++;
  135.          if (l2<=l1)
  136.          {
  137.         REPORT
  138.             i = l2-f1; while (i--) *el++ = *el1++ + *el2++;
  139.             i = l1-l2; while (i--) *el++ = *el1++;
  140.             i = l-l1;  while (i--) *el++ = 0.0;
  141.          }
  142.          else
  143.          {
  144.         REPORT
  145.             i = l1-f1; while (i--) *el++ = *el1++ + *el2++;
  146.             i = l2-l1; while (i--) *el++ = *el2++;
  147.             i = l-l2;  while (i--) *el++ = 0.0;
  148.          }
  149.       }
  150.    }
  151. }
  152.  
  153. void MatrixRowCol::Sub(const MatrixRowCol& mrc1, const MatrixRowCol& mrc2)
  154. {
  155.    int f = skip; int l = skip + storage;
  156.    int f1 = mrc1.skip; int l1 = f1 + mrc1.storage;
  157.    if (f1<f) f1=f; if (l1>l) l1=l;
  158.    int f2 = mrc2.skip; int l2 = f2 + mrc2.storage;
  159.    if (f2<f) f2=f; if (l2>l) l2=l;
  160.    real* el = store + f;
  161.    real* el1 = mrc1.store+f1; real* el2 = mrc2.store+f2;
  162.    if (f1<f2)
  163.    {
  164.       register int i = f1-f; while (i--) *el++ = 0.0;
  165.       if (l1<=f2)                              // disjoint
  166.       {
  167.      REPORT                                // not accessed
  168.          i = l1-f1;     while (i--) *el++ = *el1++;
  169.          i = f2-l1;     while (i--) *el++ = 0.0;
  170.          i = l2-f2;     while (i--) *el++ = - *el2++;
  171.          i = l-l2;      while (i--) *el++ = 0.0;
  172.       }
  173.       else
  174.       {
  175.          i = f2-f1;    while (i--) *el++ = *el1++;
  176.          if (l1<=l2)
  177.          {
  178.         REPORT
  179.             i = l1-f2; while (i--) *el++ = *el1++ - *el2++;
  180.             i = l2-l1; while (i--) *el++ = - *el2++;
  181.             i = l-l2;  while (i--) *el++ = 0.0;
  182.          }
  183.          else
  184.          {
  185.         REPORT
  186.             i = l2-f2; while (i--) *el++ = *el1++ - *el2++;
  187.             i = l1-l2; while (i--) *el++ = *el1++;
  188.             i = l-l1;  while (i--) *el++ = 0.0;
  189.          }
  190.       }
  191.    }
  192.    else
  193.    {
  194.       register int i = f2-f; while (i--) *el++ = 0.0;
  195.       if (l2<=f1)                              // disjoint
  196.       {
  197.      REPORT                                // not accessed
  198.          i = l2-f2;     while (i--) *el++ = - *el2++;
  199.          i = f1-l2;     while (i--) *el++ = 0.0;
  200.          i = l1-f1;     while (i--) *el++ = *el1++;
  201.          i = l-l1;      while (i--) *el++ = 0.0;
  202.       }
  203.       else
  204.       {
  205.          i = f1-f2;    while (i--) *el++ = - *el2++;
  206.          if (l2<=l1)
  207.          {
  208.         REPORT
  209.             i = l2-f1; while (i--) *el++ = *el1++ - *el2++;
  210.             i = l1-l2; while (i--) *el++ = *el1++;
  211.             i = l-l1;  while (i--) *el++ = 0.0;
  212.          }
  213.          else
  214.          {
  215.             REPORT
  216.             i = l1-f1; while (i--) *el++ = *el1++ - *el2++;
  217.             i = l2-l1; while (i--) *el++ = - *el2++;
  218.             i = l-l2;  while (i--) *el++ = 0.0;
  219.          }
  220.       }
  221.    }
  222. }
  223.  
  224.  
  225. void MatrixRowCol::Add(const MatrixRowCol& mrc1, real x)
  226. {
  227.    REPORT
  228.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  229.    if (f < skip) { f = skip; if (l < f) l = f; }
  230.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  231.  
  232.    real* elx = store+skip; real* ely = mrc1.store+f;
  233.  
  234.    int l1 = f-skip;  while (l1--) *elx++ = x;
  235.        l1 = l-f;     while (l1--) *elx++ = *ely++ + x;
  236.        lx -= l;      while (lx--) *elx++ = x;
  237. }
  238.  
  239. void MatrixRowCol::RevSub(const MatrixRowCol& mrc1)
  240. {
  241.    REPORT
  242.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  243.    if (f < skip) { f = skip; if (l < f) l = f; }
  244.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  245.  
  246.    real* elx = store+skip; real* ely = mrc1.store+f;
  247.  
  248.    int l1 = f-skip;  while (l1--) { *elx = - *elx; elx++; }
  249.        l1 = l-f;     while (l1--) { *elx = *ely++ - *elx; elx++; }
  250.        lx -= l;      while (lx--) { *elx = - *elx; elx++; }
  251. }
  252.  
  253. void MatrixRowCol::Copy(const MatrixRowCol& mrc1)
  254. {
  255.    REPORT
  256.    int f = mrc1.skip; int l = f + mrc1.storage; int lx = skip + storage;
  257.    if (f < skip) { f = skip; if (l < f) l = f; }
  258.    if (l > lx) { l = lx; if (f > lx) f = lx; }
  259.  
  260.    real* elx = store+skip; real* ely = mrc1.store+f;
  261.  
  262.    int l1 = f-skip;  while (l1--) *elx++ = 0.0;
  263.        l1 = l-f;     while (l1--) *elx++ = *ely++